A  Probabilistic  Multiscale  Approach  to  Hysteresis  in 
Shear  Wave  Propagation  in  Biotissue 


H.T.  Banks 

Center  for  Research  in  Scientific  Computation 
North  Carolina  State  University 
Raleigh,  NC  27695-8205 


Gabriella  A.  Pinter 
Department  of  Mathematical  Sciences 
University  of  Wisconsin,  Milwaukee 
Milwaukee,  WI  53201-0413 

January  31,  2004 


Abstract 

Motivated  by  a  problem  involving  wave  propagation  through  viscoelastic  biotissue, 
we  present  a  theoretical  framework  for  treating  hysteresis  as  multiscale  phenomena 
which  must  be  averaged  across  distributions  of  internal  variables.  The  resulting  sys¬ 
tems  entail  probability  measure  dependent  partial  differential  equations  for  which  we 
establish  well-posedness  in  a  framework  that  leads  readily  to  computationally  useful 
approximations. 

Keywords:  multiscale  modeling,  hysteresis,  wave  propagation,  biotissues,  probability  measures, 
Prohorov  metric 


1  Introduction 

At  the  heart  of  modern  day  sciences  (especially  the  computational  sciences)  of  materials, 
biology,  etc.,  is  the  challenge  of  effectively  using  the  explosion  of  available  information  at 
micro  (genomic,  atomic,  molecular,  nano,  etc.)  levels  to  develop  more  and  more  accurate 
biologically  and  physically  based  system  response  (macro  level)  models  for  use  in  prediction, 
control,  design,  and  simple  organism/structure  simulations.  A  number  of  fundamental  issues 
(e.g.,  how  to  model  uncertainty /variability  in  heterogeneous  materials)  in  multiscale  model¬ 
ing  and  control  are  important  in  addressing  this  challenge.  In  this  presentation  we  discuss 


1 


Report  Documentation  Page 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1.  REPORT  DATE 

31  JAN  2004  2.  REPORT  TYPE 

3.  DATES  COVERED 

00-00-2004  to  00-00-2004 

4.  TITLE  AND  SUBTITLE 

A  Probabilistic  Multiscale  Approach  to  Hysteresis  in  Shear  Wave 
Propagation  in  Biotissue 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

North  Carolina  State  University, Center  for  Research  in  Scientific 
Computation,  Raleigh,  NC, 27695-8205 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

see  report 

15.  SUBJECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF:  17.  LIMITATION  OF 

18.  NUMBER  19a.  NAME  OF 

a.  REPORT  b.  ABSTRACT  c.  THIS  PAGE 

unclassified  unclassified  unclassified 

20 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


some  of  these  issues  in  the  context  of  internal  dynamics  (molecular  level)  in  system  response 
models  (organism  level).  More  precisely,  we  offer  a  probabilistic  based  internal  strain  variable 
framework  to  treat  nonlinear  hysteresis  arising  in  propagation  of  shear  waves  in  biotissues. 
The  usual  models  for  such  problems  involve  integro-partial  differential  equations  which  are 
most  often  phenomenological  in  nature  as  well  as  being  computationally  challenging.  Our 
approach  is  founded  on  the  belief  that  hysteresis  is  actually  a  manifestation  of  the  presence 
of  multiple  scales  in  a  physical  or  biological  system  that  is  frequently  modeled  (and  masked) 
with  a  phenomenological  representation  such  as  an  hysteresis  integral  for  the  macroscopic 
stress-strain  constitutive  law  (  e.g.,  for  models  of  viscoelastic  damping  via  hysteretic  integrals 
versus  internal  variable  representations,  see  [13]  and  the  references  therein).  The  approach 
we  present  embodies  a  multiscale  treatment  of  hysteresis  as  resulting  from  a  continuum  of 
microscopic  or  molecular  level  internal  strain  mechanisms  averaged  (with  respect  to  a  mate¬ 
rial  dependent  probability  distribution)  across  populations  of  molecules.  Such  an  approach 
using  nonlinear  stick-slip  molecular  reptation  models  for  both  tensile  and  shear  deforma¬ 
tions  has  recently  [11], [12]  been  developed  and  successfully  used  in  computational  efforts 
with  composite  rubbers.  The  probabilistic  framework  we  present  here  has  its  origins  in  work 
on  distributed  growth  rates  for  marine  (mosquitofish)  populations  in  size-structured  models 
as  developed  in  [1] ,  [5] ,  [6] ,  [7] ,  [8] .  We  develop  here  for  the  first  time  a  theoretical  framework 
for  molecular  based  strain  models  which  readily  leads  to  efficient  computational  methods 
that  offer  an  alternative  to  the  computationally  intensive  phenomenological  viscoelastic  soft 
tissue  kernels  proposed  by  Fung  and  his  colleagues.  In  [22] ,  Fung  lays  the  foundations  for  the 
essentials  of  polymer  biomaterials  such  as  clastin,  fibers,  collagen,  etc.,  which  experimentally 
manifest  classical  hysteretic  behavior  in  tension  and  shear  (see  p.  261,  [22]).  He  furthermore 
cogently  argues  the  essential  nonlinear  nature  of  these  materials.  He  also  explains  the  need 
for  continuous  spectrum  models  (e.g.,  a  continuum  of  relaxation  parameters  )  due  to  the 
strain  rate  independence  of  the  stress-strain  law  ([22]  p.  281-287)  over  several  decades  of 
frequencies  (see  also  [21],  [25],  [26],  [28]).  This  led  to  the  quasi-linear  Fung  kernel  (  see  (2.6) 
below)  in  the  context  of  Boltzman  type  model  (2.5)  constitutive  laws.  The  formulation  we 
give  here  also  provides  continuous  spectrum  models,  although  we  subsequently  report  on 
calculations  in  [2]  with  finite  spectrum  approximate  models  (Section  6  below).  The  com¬ 
putational  results  suggest  that  such  models  will  often  suffice  to  describe  well  the  type  of 
deformations  arising  in  the  wave  propagation  problems  motivating  our  developments  in  this 
paper. 

Our  motivating  problem  arose  in  joint  collaborations  with  scientists  and  engineers  at 
MedAcoustics,  Inc.,  as  part  of  the  Industrial  Applied  Mathematics  Program  at  N.C.  State 
University  (  see  http://www.ncsu.edu/crsc/  ).  We  summarize  briefly  here  the  salient  points 
of  this  problem,  while  referring  the  reader  to  [2]  and  [14]  for  more  details.  Turbulent  blood 
flow  due  to  arterial  stenosis  in  partially  occluded  arteries  produces  normal  forces  on  arterial 
walls.  This  results  in  vibrations  in  the  surrounding  body  tissue  which  are  transmitted  to 
body  surfaces  in  two  forms:  compressional  waves  and  shear  waves.  The  shear  waves  are  at 
low  frequencies  (  <  2kHz  )  with  low  propagation  speed  and  attenuation.  Devices  involving 
multiple  arrays  of  piezoceramic  sensors  were  developed  at  MedAcoustics,  with  the  goal  of 
measuring  shear  wave  propagation  at  the  surface  of  the  chest.  The  resulting  signals  are  then 
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processed  to  determine  the  location  and  extent  of  the  stenosis.  A  part  of  this  overall  detection 
and  diagnostic  problem  is  the  focus  of  our  efforts:  modeling  of  the  shear  wave  propagation 
through  a  heterogeneous  chest  cavity.  The  cavity  is  composed  of  soft  body  tissues  (lung, 
muscular  and  connective  tissue)  as  well  as  bone.  Here  we  consider  shear  wave  propagation 
through  a  viscoelastic  heterogeneous  medium.  At  MedAcoustics,  early  experiments  and  data 
collection  were  carried  out  on  cylindrical  shaped  mold  of  tissue-like  synthetic  gel  as  depicted 
schematically  in  Figure  1.  The  cylinder  of  gel  surrounded  a  tube  in  which  source  disturbances 
simulating  disrupted  flows  were  produced.  Thus  the  tube  with  disturbances  mimicked  an 
artery  with  stenosis.  Shear  waves  were  measured  with  sensor  arrays  mounted  on  the  gel  outer 
surface.  Subsequent  experiments  involved  data  collection  on  pigs  (who,  unfortunately  for 
them,  have  cylindrical-like  chest  geometries!!).  In  [2],  the  authors  developed  and  used  a  one¬ 
dimensional  (axial  symmetry  was  assumed)  homogeneous  medium  model  for  computations 
and  analysis.  Here  we  consider  an  extension  of  that  model  which  allows  for  “molecular”  level 
heterogeneity  in  the  tissue  (gel),  while  retaining  the  overall  geometry  and  axial  symmetry. 
An  inner  radius  R\  and  outer  radius  R2  for  the  gel  is  assumed,  with  the  gel  initially  at  rest. 
The  outer  surface  of  the  gel  is  a  free  surface. 


Acoustic  disturbance 


Figure  1:  Cylindrical  model  geometry  with  artificial  “stenosis”  and  piezo  sensors 


2  Problem  formulation 

In  [2]  and  [14]  the  following  one- dimensional  model  (see  Figure  1)  is  introduced  to  describe 
the  propagation  of  shear  waves  in  soft  tissue 

d2u  d 

p-^  -  g^a(t,x)  =  F(t,x),  Ri  <  x  <  R2, 
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where  p  is  the  mass  density,  u  is  the  Cartesian  shear  displacement,  a  is  the  shear  stress  and 
F  represents  a  body  forcing  term.  We  assume  that  at  the  inner  or  left  boundary,  (x  =  R\), 
we  have  a  pure  shear  shearing  force,  while  the  outer  or  right  boundary,  (x  =  R2),  is  a  free 
surface 


cr(t,Ri)  =  fit),  and  a(t,R2)  =  0. 


The  initial  conditions  are 

w(0,x)  =  uo(x),  and  ut(0,x)  =  u\(x),  R\  <  x  <  R2. 

To  complete  the  model  one  needs  to  find  an  adequate  constitutive  relationship  that  relates 
the  shear  stress,  a,  and  strain,  £  =  ux(t,x ),  in  soft  tissue  (arteries,  muscle,  skin,  lung,  etc., 
. . . ).  In  [2]  an  internal  strain  variable  model  is  considered  in  the  following  form 

N 

a(t,x)  =  CDi(t,x)  +  ^2£j(t,x),  (2.1) 

3= 1 

with 

ds At )  1  .  .  d  .  ,  w 

——  =  -- £j(t )  +  Cj-ae{ux{t)) 

£j(0,x)  =  0,  j  =  l,...N, 

where  oe  is  the  elastic  response  function  defined  in  Fung  ([22],  §7)  and  given  as 

ae{ux)  =  !  +  (2.4) 

On  the  other  hand,  Fung  proposes  the  constitutive  relationship 

ait)  =  [  Git  —  s)^j— (2.5) 
Jo  ds 

where  G(t)  is  the  reduced  relaxation  function  given  in  the  form 

G(t)  =  |l  +  C[Ei( — )  —  Ei(  — )]}  [1  +  cln(  — )]  1.  (2-6) 

l  t2  tx  J  tx 

Here  E\(z)  =  ff°  ^ dt ,  C  represents  the  degree  to  which  viscous  effects  are  present,  and  T\ 
and  t2  represent  the  fast  and  slow  viscous  time  phenomena.  We  note  that  the  internal  strain 
variable  formulation  (2.1)-(2.3)  is  equivalent  to  the  constitutive  relationship  proposed  by 
Fung  if  one  considers  an  approximation  of  the  relaxation  function  G  by  a  sum  of  exponential 
terms.  Various  internal  strain  variable  models  are  investigated  in  [2]  and  a  good  agreement 
is  demonstrated  between  the  two  internal  strain  variable  model  (cr  =  eq  -\-£2)  and  undamped 
simulated  data  based  on  the  Fung  kernel  G.  We  also  remark  that  theoretical  well-posedness 
results  for  the  one  strain  variable  model  (which  can  be  readily  extended  to  treat  finite  multiple 


(2.2) 

(2.3) 
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strain  variables)  are  given  in  [14].  We  further  note  that  the  derivations  in  [11],  [12],  and  the 
references  therein  yield  a  molecular  based  foundation  for  internal  variable  models  such  as 
(2.2)-(2.3).  These  are  generalizations  of  the  linear  stick-slip  models  of  Doi  and  Edwards  [20] 
and  Johnson  and  Stacer  [23]. 

As  explained  in  the  Introduction,  here  we  are  concerned  with  a  generalization  of  the 
internal  strain  variable  formulation.  Fung’s  proposed  model  provides  a  continuous  spectrum 
of  relaxation  times  in  contrast  with  the  finite  multiple  internal  strain  model.  In  the  latter, 
the  material  is  assumed  to  have  discrete  relaxation  rate  constants,  Tj,  that  correspond  to 
a  discrete  hysteresis  spectrum.  As  we  have  noted,  Fung  argues  that  this  is  incompatible 
with  the  observation  that  the  hysteresis  loop  is  independent  of  the  strain  rate  within  several 
decades  of  the  rate  variation  for  many  biological  soft  tissues.  Here  we  propose  to  realize  a 
continuous  spectrum  of  relaxation  times  in  the  form 

a(t,x;P)  =  CDe(t)  +  f  G(t  -  s,P)^-ae{ux(s))ds,  (2.7) 

Jo  as 

with 

G(t,P)  =  Jrg(t,r)dP(T),  (2.8) 

where  T  —  [ri,  72]  C  (0,  00),  P  is  a  probability  density  function  on  T,  and  g  is  a  continuous 
function  of  relaxation  times  r  on  T.  We  note  that  the  case  g(t,r)  =  e ~rt  corresponds  to  a 
continuum  of  internal  strain  variable  models  “weighted”  by  the  probability  density  function 
P.  This  can  be  seen  readily  by  applying  the  variation  of  constants  formula  in  (2.2)-(2.3) 
(with  Cj  =  1,  j  =  1, . . . ,  N) 

,,  \  f*  t-s )  d  ,  /  w  1 

£j{t',Ti)=  /  e  Tj  —cre(ux(s))ds, 

Jo  as 

<j(t,x;P)  =  CDux{t)  +  [  [  e  T>(t  A)  ^-ae{ux(s))ds  dP(r), 

Jt  Jo  ds 

(and  this  can  be  generalized  using  a  product  measure  on  T  x  C  to  include  a  continuum  of 
Cj- s).  In  the  next  section  we  present  existence-uniqueness  results  for  the  generalized  model 
with  constitutive  relationship  (2.7)-(2.8). 

3  Existence  and  uniqueness  of  weak  solutions 

Motivated  by  the  weak  or  variational  form  of  the  model  derived  in  the  previous  section  we 
consider  the  abstract  system 


ua  -  x\P)  =  F  in  V* 

(3.1) 

u( 0,  x)  =  u0  G  V  =  i71(i?i,  R2) 

(3.2) 

ut( 0,x)  —  e  H  —  L2(Ri,R2) 

(3.3) 
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with  boundary  conditions 


a(t,RuP)  =  f(t) 
nit .  Rz,  P)  =  0, 


(3.4) 

(3.5) 


where  the  shear  stress  a  is  given  as 

ft 

a(t,x]P)  —  CDux(t,x)  +  /  G(t- s,P)—cre(ux(s,x))ds,  (3.6) 

Jo  CIS 

with 

G(t,P)  =  Jrg(t,T)dP(r)>  (3.7) 

and  P  is  a  probability  density  function  on  T.  We  will  denote  the  norm  and  inner  product  in  H 
by  ||  •  ||  and  (•,  •),  while  all  other  norms  and  inner  products  will  be  appropriately  labeled.  We 
note  that  \\v\\v  =  IMI2  +  ||w||2,  and  that  V  is  compactly  embedded  in  H  with  V  H  V* 
forming  a  Gelfand  triple  [29]. 

We  make  the  following  assumptions 

(Al)  The  forcing  term  F  satisfies  F  E  L2( 0,  T;  V*). 

(A2)  The  inner  boundary  condition  satisfies  /  G  L2(0,T). 

(A3)  The  elastic  response  function  ae  satisfies  the  following  local  Lipschitz  condition 

\\ae(u)  -  ae{v)\\  <  LBr\\u  -  v\\ 

for  some  positive  constant  L#r  and  for  all  u,v  E  Bh( 0,  r),  the  ball  in  F[  centered  at  0 
of  radius  r. 

(A4)  There  exist  positive  constants  C\  and  C-2  such  that 

||cre(w)||  <  Ci|H|  +  C2 


for  every  u  E  H. 

We  remark  that  these  are  the  same  assumptions  as  those  of  [14].  In  addition,  we  make 
the  following  assumption  on  the  relaxation  function  G(t,  P )  =  Ir9(t  -  s,  r)dP(r). 

(A5)  The  function  g  is  continuous  in  r  on  T  =  [ti,t2]  C  (0,  oo),  g  is  differentiable  with 
respect  to  t  on  ]R+ ,  and  there  exist  constants  C3  and  C4  such  that  [<?(£,  t)|  <  G3  and 
| g(t,  r)|  <  C\  for  all  t  E  1R+  and  all  r  E  T. 

We  note  that  g(t,r)  =  e ~^t,  which  corresponds  to  a  continuum  of  internal  strain  variable 
models,  satisfies  this  latter  condition. 

We  define  weak  solutions  for  the  system  (3.1)-(3.7)  in  the  following  way. 
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Definition  3.1  Let  Ct  =  {w  :  [0,T]  — >  H,  w  G  CV(0,T;fo)  fl  L2(0,T]V)  and  wt  G 
Cw{ 0,  T;  H)  fl  L2( 0,  T;  fo)}.  We  say  that  u  &  Ct  is  a  weak  solution  of  the  system  (3.1)-(3.7) 
if  it  satisfies 

rt  rs  (l 

/  [-(us(s),T]s(s))  +  CD(usx(s),r]x(s))  +  (  /  G(s  -  u,  P)—ae(ux(s))dv,  r)x(s))]ds 
Jo  Jo  as 

+  (ut(t),rj(t))  -  <i/i, 77(0))  =  [  [( F(s),rj(s))v*y  -  f(s)rj(s,Ri)]ds 

Jo 

for  any  t  G  [0,  T]  and  77  G  Ct  with  initial  conditions  uq  G  V  and  ru,\  G  H. 

As  in  [14]  we  note  that  this  notion  of  the  weak  solution  for  system  (3.1)-(3.7)  agrees  with 
the  usual  one  in  that  it  yields  utt  G  L2(0,T;  V*)  with  equation  (3.1)  holding  in  the  sense  of 
L2(0,T;  V*). 

Our  first  result  is  that  the  system  (3.1)-(3.7)  has  a  unique  weak  solution. 


Theorem  3.1  Under  assumptions  (Al)-(A5)  the  system  (3.1)-(3.7)  has  a  unique  global 
weak  solution  on  any  finite  interval  [0,  T], 

Proof:  Our  arguments  to  establish  this  statement  essentially  follow  those  of  [14].  Thus  we 
present  only  a  short  outline  and  point  out  the  similarities  between  the  one  internal  variable 
model  of  [14]  and  its  generalization  considered  here.  First  we  note  that  the  constitutive 
relationship  a  =  Cr>ux  +  E\  that  was  used  in  [14]  leads  to 


cr(t)  =  CDux(t)  +  J  e  s)—  ae(ux(s))ds 

=  CDiix(t)  +  ae{ux{t))  -  e~Aae(u0x)  -  /  —  (e~^{t~s))ae(ux(s))ds, 

Jo  ds 


while  here  we  propose 


rt  (J 

cr(t)  =  CDux{t)  +  /  G(t-  s,P)—ae(ux(s))ds 
Jo  ds 

=  CDux(t)  +  G(0,P)ae(ux(t))  -  G(t,  P)ae(u0x)  +  f  G(t  -  s,P)ae(ux(s))ds, 

Jo 

with 

G(t,P)  =  J^g(t,T)dP(T). 

Thus  as  long  as  G(t,  P)  has  the  properties  that  were  used  in  [14]  in  connection  with  the 
exponential  kernel,  we  can  essentially  repeat  all  the  arguments  there  and  can  arrive  at  the 
desired  conclusion.  We  remark  that  assumption  (A5)  guarantees  this  analogy  between  the 
two  problems.  Hence  we  can  proceed  as  in  [14]  and  as  a  first  step  establish  the  local  existence 
of  weak  solutions  under  the  assumptions  (Al)-(A3)  and  (A5).  To  this  end  we  define  the 
radial  retraction  Pi  from  H  onto  Bh(uqx,  1)  and  let  cre  =  ae(P\U ),  for  u  G  H.  We  note  that 
cre  satisfies  a  global  Lipschitz  property 


ae(u)  —  (Te(u)||  <  L\\u  —  v||  for  all  u,v  G  p[, 
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and  a  boundedness  property 


||de(w)||  <  Ci|M|  +  C2 

for  all  u  G  H  with  positive  constants  C\  and  C2. 

Next,  we  let  {^}“1  be  any  linearly  independent  total  subset  of  V,  and  for  any  m,  m  = 
1, . . we  form  Vrn  =  span{^ i, . . .  We  choose  sequences  {w™}  and  {u™}  with  u™  in 

Vm,  i  =  0,1,  such  that  n™  — >  u0  in  V  and  n™  — >  Ui  in  H  as  m  — >  oo.  We  define  the  Galerkin 
approximation 

m 

fc=i 

as  the  unique  solution  of  the  m-dimensional  integro-differential  system 

+  CD(u^x,^jx)  +  G(t  —  s,  P)-^ae(u™(s))ds,  il)jx) 

=  (F(t),ipj)v*,v  ~  f(t)il>j(Ri)  (3.8) 

for  j  =  1 , ...  ,m  on  the  interval  [0,  T],  for  some  T  >  0.  Now  we  can  use  the  same  techniques 
as  in  [14]  to  establish  an  a  priori  estimate 

\Km2  +  \Km2  +  CD  f\K(s)fdS<K,  t  e  [0,T],  (3.9) 

Jo 

where  K  =  K{uq,U\,  f,  F,Cd,T,  &e,  g,  P)  and  is  independent  of  m.  Thus  we  can  conclude 
that  there  exists  a  function  u  G  L2(0, T;  V)  such  that 

um  — "  u  weakly  in  L2(0,  T;  V) 


and 

■u™  — ^  ut  weakly  in  L2(0,  T;  V). 

The  following  additional  convergences  that  play  a  crucial  role  in  the  proof  can  also  be  es¬ 
tablished  by  the  same  methods  as  those  used  in  [10,  14]  (Ascoli-Arzela  Theorem,  Aubin’s 
lemma) . 

(Rl)  um{t)  — ^  u(t)  weakly  in  V  uniformly  in  t  G  [0,  T\,  i.e.,  um  — ■>  u  in  Cw(0,  T;  V)). 

(R2)  u™{t)  — ^  ut(t)  weakly  in  F[  uniformly  in  t  G  [0,  T],  i.e.,  w™  — >  ut  in  CV(0,  T;  H)). 

(R3)  u™  — »  ut  in  L2(0,  T;  H). 

(R4)  There  exists  a  function  h  G  L2(0,  T;  H)  such  that 

/  G(t  —  s)aJux(s))ds  — *■  h 
Jo 


weakly  in  L2(0,T;  H). 


We  remark  that  these  convergences  also  hold  on  any  subinterval  of  [0,  T\.  Now  we  verify  that 
u  is  a  weak  solution  of  (3.1)-(3.7)  with  ae  instead  of  ae.  First  we  define  the  class  of  functions 
Bm  =  {q  G  FT |  rj(t)  =  J2kLick(t)i/>ki  ck  £  C,1[0,T]},  and  note  that  B  =  U m=1Bm  is  dense 
in  CT.  Now  we  multiply  (3.8)  by  Ck(t),  sum  from  1  to  M  and  integrate  over  (0,  t)  to  obtain 
that  for  all  rj  G  Bm,  m>  M 

f  [«s(s),v(s))v*y  +  CD(u™{s),qx{s ))  +  ([  G(s-u,  P)ae(u™(u))du,  qx{s))}ds 
Jo  Jo 

=  f  {F(s),q(s))v*y  -  f(s)q(s,R1)ds. 

Jo 

We  integrate  the  first  term  by  parts  and  then  take  the  limit  as  m  — >  oo  using  the  convergences 
(R1)-(R4)  to  arrive  at 


[  [-(us(s),qs(s))  +  CD(usx(s),qx(s)}  +  (h(s) ,  rjx(s))]ds 
Jo 

+(ut(t),ri(t))  -  (uyrjiO)) 

=  [  (F(s),ri(s))vy  -  f(s)v(s,Ri)ds 
Jo 

for  all  q  G  Ft-  Thus  we  can  claim  that  u  is  a  weak  solution  of  (3.1)-(3.7)  (with  ae  replaced 
by  <xe)  if  we  show  that 

rt  rt  rs 

/  (h(s),qx(s))ds  =  /  (/  G(s  —  is,  P)de(ux(is))dv,  qx(s))ds 
Jo  Jo  Jo 

for  all  q  G  Ft-  This  is  achieved  by  establishing  the  strong  convergence  u™(t)  — >  ux(t )  in  H  as 
m  — »  oo,  again  using  the  same  arguments  as  in  [14].  Uniqueness  of  weak  solutions  is  proved 
the  standard  way. 

Now  we  claim  that  the  original  system  (3.1)-(3.7)  with  cre  has  a  unique  local  weak  solution. 
Since  ux  is  continuous  in  t  ([19],  p.555)  it  follows  that  there  exists  t* ,  0  <  t*  <  T  such  that 

||'lta;(f)  WQa;||  ^  1 


for  all  t  G  [0,  £*].  Therefore 


Geiuxit))  =  ae(ux(t)),  t  G  [0 ,t*], 

i.e.,  u  is  a  weak  solution  of  (3.1)-(3.7)  on  [0,  t*]. 

Finally,  we  can  conclude  the  proof  of  Theorem  3.1  by  showing  that  this  local  solution  actu¬ 
ally  exists  on  any  arbitrary  interval  [0,  T],  This  is  accomplished  by  utilizing  the  boundedness 
assumption  (A4)  to  obtain  an  a  priori  estimate  similar  to  (3.9)  for  Galerkin  approximations 
involving  ae  instead  of  ae,  and  then  argue  that  ux(t)  is  pointwise  bounded.  The  local  Lips- 
chitz  condition  (A3)  can  then  be  applied  and  similar  steps  to  those  above  can  be  repeated 
to  show  that  u  is  the  unique  weak  solution  of  (3.1)- (3. 7). 
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4  Continuous  dependence 

In  this  section  we  show  that  the  unique  weak  solution  of  system  (3.1)-(3.7)  depends  continu¬ 
ously  on  the  “distributed  relaxation  times” .  This  result  is  important  for  our  ultimate  goal  of 
identifying  stenosis,  i.e.,  in  considerations  involving  the  inverse  problem.  To  this  end,  we  first 
introduce  a  topology  on  the  probability  distribution  functions  (PDFs)  P  G  T2]),  where 

7,([ti,t2])  is  the  set  of  all  PDFs  on  the  Borel  subsets  of  A  convenient  topology  that 

will  also  be  useful  in  treating  the  general  inverse  problem  is  given  by  the  Prohorov  metric 
[3,  4,  15].  For  the  sake  of  completeness  we  give  the  basic  definitions  and  results  concerning 
this  metric  here. 

Suppose  that  ( Q,d )  is  a  complete  metric  space,  and  for  any  closed  subset  F  C  Q  and 
e  >  0  define  the  set 

Fe  =  {qeQ  \  d(q,  q)  <  e,  g  G  F}. 

Then  the  Prohorov  metric  p  :  V(Q)  x  V(Q)  —■ ►  JR+  is  defined  by 

p(Pi,  P2)  =  inf{e  >  0  :  P^F]  <  P2[Fe }  +  e,  F  closed,  F  C  Q}. 

It  is  well-known  that 

(1)  ( V(Q),p )  is  a  complete  metric  space,  and 

(2)  if  Q  is  compact  then  ( V(Q),p )  is  a  compact  metric  space. 

The  crucial  properties  of  this  metric  in  connection  with  our  continuous  dependence  and 
inverse  problems  are  summarized  in  the  following  theorem  ([15]). 

Theorem  4.1  Given  Pk,  P  G  V(Q),  the  following  convergence  statements  are  equivalent 

(a)  p(Pk,P )  0; 

(b)  fQ  fdPk(q)  —>  fQ  fdP(q )  for  all  bounded,  uniformly  continuous  functions  f  :  Q  — >  M; 

(c)  Pk[A]  — >  P[A]  for  all  Borel  sets  A  C  Q  with  P[dA]  =  0; 

We  can  now  proceed  to  show  that  the  weak  solution  of  system  (3.1)-(3.7)  depends  con¬ 
tinuously  on  the  probability  distribution  function  P  in  the  Prohorov  metric  topology. 

Theorem  4.2  Assume  that  assumptions  (A1)-(A5)  hold.  If  Pk  — >  P  in  the  Prohorov  metric 
p,  then  Uk  — >  u  in  L2( 0,  T;  V)  and  Ukt  — 1 >  ut  in  L2( 0,  T;  V),  where  Uk  are  the  weak  solutions 
corresponding  to  probability  density  functions  Pk  and  u  is  the  weak  solution  corresponding  to 
P. 

Proof:  We  first  note  that  if  Pk  —■ ►  P  and  G  is  defined  as  G(t,  P )  —  It  g(t,r)dP(r ),  where  g 
satisfies  (A5),  then  by  Theorem  4.1.  we  have  that 

G(t,  Pk)  — > >  G(t,  P )  for  every  t  >  0, 

G(t,  Pk)  — > >  G(t,  P)  for  every  t  >  0. 
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(4.1) 

(4.2) 


We  also  remark  that  (A5)  implies  that 


(4.3) 

(4.4) 


\G(t,Pk)\  <  C:h  for  every  k,  k  =  1, . . . 

\G(t,Pk)\  <  C4,  for  every  k,  k  —  1, . . . 

Now  we  have  that 

a(t,  x;  P)  =  CDe  +  f  G{t  -  s)^-(cre{ux{s)))ds  =  CDi  +  Ei(t,£]P), 

Jo  as 

where  e  =  ux,  and  by  the  above  remarks  Ei(t,  e\  Pk)  — >  Ei(t,  £,  P )  for  a  hxed  s,  and  any  t  >  0 
whenever  Pk  — >  P  in  the  Prohorov  metric.  Let  u™ ,  um  denote  the  Galerkin  approximations 
corresponding  to  the  weak  solutions  Uk  and  u  of  (3.1)-(3.7)  with  respective  PDFs  Pk  and  P. 
We  note  that  the  initial  conditions  do  not  depend  on  the  PDFs  Pk  or  P.  Thus  we  have 

{uZMv-y  +  CD{u^jx)  +  <Pi«;  P)^jx)  =  (F(t),^)v,y  -  ,) 


and 


(uktv  'lPj)v*,v 


+  CD(u 


m 

hix") 


which  yields 


Ku  -  <C <Pj )v-,v  +  CD(u%x  -  <,  fe)  +  P)  -  P),lW  =  0-  (4.5) 

Let  AE™  =  Ei(u™x]  P)  —  Ei(u™’,  P).  We  remark  that  by  the  above  observation  on  the 
boundedness  of  G  and  G  we  have  that  the  same  a  priori  estimates  can  be  shown  to  be 
valid  for  both  the  Galerkin  approximations  {u™}  and  {um},  i.e.,  there  exists  a  constant  K 
independent  of  m  and  k  such  that 

\KM\2  +  \KM\2  +  cD  f  \Ksx(s)\\2ds  <  k 

\K(t)\\2  +  \K(t)\\2  +  CD  f  \\uZ(s)\\2ds  <  K, 

Jo 

for  all  t  G  [0,T],  m  —  1, . . . .  Since  these  estimates  guarantee  the  pointwise  boundedness  of 
\\u™\\  and  1 1 ux |  (by  the  weak  lower  semicontinuity  of  norms  and  the  convergence  u™(t )  — > 
ux(t)  in  H),  we  can  utilize  the  local  Lipschitz  property  (A3)  of  ae  in  Bh(0,  VK)  together 
with  the  growth  condition  (A4)  to  obtain  the  following  estimate  for  AE™ 


\\AE™it)\\  <  \G(0,Pk)\L^\\uZ(t)  ~u™(t)\\  +  |G(0,Pfc)  -  G(0,  P)|{G1|K(t)||  +  G2} 
+  \G(t,Pk)\L^\\u^M-u™m  +  \G(t,Pk)  ~G(f,P)|{C1|K(0)||  +C2} 

+  J‘{\G(t-s,Pt)\Lv]!\KM-u'Z(S)\\ 

+|G(t  -  s,  P„)  -  CAt  -  s,  P)|{C,  ||«s)  +  C2}]*. 
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Let  A ™(t)  =  u™(t)  —  um(t).  The  a  priori  estimate  yields  that  there  exist  positive  constants 
Mi,  M2  and  M3  such  that 

l|ASr(«)lt  <  Mi||A£(<)||  +  M,\G(0,Pt)  -  G(0,P)\  +  M3\G(t,Pk)  -  G(t,P) | 

/'[M2 ||  AJ”  (»)  ||  +  M3\G(t  -  s,  Pk)  -  G(t  -  s,  P)\]ds,  (4.6) 

Jo 

where  we  used  the  fact  that  u™x( 0)  —  u™( 0)  =  0.  We  now  return  to  (4.5),  multiply  by 
-  Ojl(t)  and  sum  over  j  to  obtain 

~||  AS(«)||2  +  CD||AS,(«)f  +  (AJ3fW.  AS„(())  =  0. 

We  add  (A™x(t),  A™tx(t))  to  both  sides,  and  then  using  standard  inequalities  we  arrive  at 

\jt  (||AS(«)H2  +  l|A£(i)f)  +  ^IIAS^t)!!2  <  ^  (l|ASr(«)H2  +  l|A£(i)H2)  ■  (4.7) 

Let  us  introduce  the  notation 

A  G{t,k)  =  G{t,Pk)-G{t,P) 

and 

A  G(t,k)  =  G(t,Pk)-G(t,P). 

By  the  estimate  (4.6)  we  have 

||ABr(*)H2  <  {||  A£(i)||2  +  |AG(0,  k) |2  +  |AG(t,*)|2  +  tj‘  ||A£(S)||2*. 

+  J‘\AG(t-s,k)\2ds\.  (4.8) 

Next  we  integrate  (4.7)  from  0  to  t  and  obtain  after  using  standard  methods  and  the  estimate 
(4.8)  that 

\\AZ(t)\\2  +  \\^TM\2  +  CD  f\\AZAs)fds<Ki  fllA^sWds  +  Mt),  (4.9) 

Jo  Jo 

where  Jkit)  is  given  by 

Jk(t)  =  ^  A|AG(0,  k) I2  +  \AG(s,  k)\2  +  %\A G(t  -  5,  k)\2}ds. 

G D  Jo  2 

By  (4.1)-(4.2)  and  (4.3)-(4.4)  we  can  assert  that 

Jkit)  0  (4.10) 

for  t  G  [0,  T]  as  k  — >  oo.  Application  of  Gronwall’s  lemma  in  (4.9)  yields 

l|AS(«)f  +  l|Al“(i)H2  <  Jkit)  +  K ,  f  Jk(s)eK'i‘-idS,  t  e  [0 ,T],  (4.11) 

Jo 


12 


By  the  weak  lower  semicontinuity  of  norms,  the  convergences  (Rl)  and  (R3)  for  a  fixed  k 
and  (4.10)  we  obtain  from  (4.11)  that 

Ukt{t)  —mt(t)  in  H  for  each  t  G  [0,  T], 


and 


We  note  that  since 


implies 


Ukxif)  ux(t)  hr  H  for  each  t  G  [0,T], 


||A™(t)-A-(0)||  <  A  ASII* 

Jo 


IiAr(t)<iiAno)n+  A  ash*, 

Jo 

which  tends  to  0  as  k  — >  oo,  we  also  have  that 

Uk(t)  — > ►  u(t)  in  H  for  each  t  G  [0,  T\. 

Additionally,  we  can  show  that  Uktx  utx  in  L2(0,T;  H)  by  using  the  estimate  (4.9).  Thus 
we  can  conclude  that  Uk  — > ►  u  in  L2(0,T;  V)  and  Ukt  — >  ut  in  L2(0,T;  V),  which  establishes 
Theorem  4.2.  We  note  that  the  proof  actually  shows  that  the  map  P  — >  u(t;  P)  is  continuous 
from  V(Q)  to  V  for  each  t  G  [0,  T], 

5  The  inverse  problem 

Since  the  final  goal  of  the  original  investigation  was  the  detection  of  stenosis  in  arteries, 
in  this  section  we  turn  to  the  study  of  the  inverse  problem.  In  particular,  we  consider 
whether  it  is  possible  to  determine  mechanical  properties  of  the  tissue  (i.e.,  the  distribution 
of  relaxation  times)  by  measuring  shear  wave  propagation  at  the  surface  of  the  chest  wall. 
Thus  we  attempt  to  estimate  P  G  V(Q),  (where  V(Q)  is  the  set  of  all  probability  density 
functions  over  Q),  given,  e.g.,  shear  displacement  data  <A,  i  —  1  . . .  ,n,  by  minimizing  the 
cost  functional 


J(P)  =  YJHU,R2-,P)-d>\2  (5.1) 

i= 1 

over  all  P  G  V(Q).  We  assume  that  Q  is  compact.  Our  continuous  dependence  result  in 
the  previous  section  (i.e.,  P  — >  u(t]P)  is  continuous  from  V(Q)  to  V,  Q  —  T,  compact) 
implies  that  the  map  P  — >  J(P)  is  continuous  from  V{Q)  to  JR.  Since  the  compactness  of  Q 
guarantees  that  of  V(Q),  we  have  that  minimizer  must  exist,  i.e.,  the  inverse  problem  has  a 
solution.  The  framework  developed  in  [3]  also  provides  a  way  to  approximate  solutions  of  the 
inverse  problem  that  is  useful  in  actual  computations.  It  is  shown  that  an  arbitrary  element 
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of  V{Q)  can  be  approximated  by  a  finite  linear  combination  of  Dirac  measures.  To  be  more 
precise,  let  Qm  —  M  =  1,2,...,  oo,  be  subsets  of  Q  chosen  so  that  Qd  =  U m=iQm 

is  dense  in  Q.  Next  define 

M 

Vm(Q)  =  {Pe  V{Q)\P  =  5>,V,  yf  e  QmiPj  rational, YPj  =  1}, 

3= i  J 

where  5q  is  the  Dirac  measure  with  atom  at  q.  Then  in  [3]  it  is  shown  that  any  element  of 
V{Q)  can  be  approximated  arbitrarily  closely  (in  the  Prohorov  metric)  by  elements  in  VM(Q) 
for  M  sufficiently  large.  That  is,  U m=iVm  (Q)  is  dense  in  V(Q)  in  the  Prohorov  metric.  In 
terms  of  our  problem  of  shear  wave  propagation  this  means  that  we  can  approximate  with  a 
model  involving  a  finite  number  of  internal  strain  variables  in  the  numerical  computations. 
The  results  of  such  an  attempt  are  detailed  in  Section  6. 

However,  in  some  problems  it  is  natural  to  minimize  J(P)  over  a  subset  of  V(Q)  consisting 
of  only  absolutely  continuous  probability  distributions,  i.e.,  those  that  have  ‘continuous’ 
densities  with  Prob[a,  b]  =  ['*  f(x)dx.  For  example,  as  we  mentioned  before,  Fung  argues 
against  a  finite  spectrum  of  relaxation  times  (which  results  if  one  uses  an  approximation 
with  finite  linear  combinations  of  Dirac  measures)  and  suggests  a  continuous  distribution  of 
relaxation  times.  In  this  context  it  might  be  more  desirable  to  minimize  J(F)  over  some 
VAQ)  —  {F1  £  F(Q)\F  =  /  /,  /  G  F},  where  F  is  a  given  set  of  functions,  and  this 
motivates  the  development  of  new  theoretical  and  approximation  results.  We  remark  that 
a  crucial  property  in  the  Banks-Bihari  [3]  framework  is  that  the  compactness  of  Q  implies 
compactness  of  the  space  ( V(Q,p ),  where  p  denotes  the  Prohorov  metric.  Thus  we  look  for 
some  type  of  compactness  condition  imposed  on  the  set  of  functions  F  that  leads  to  the 
compactness  of  Vj?(Q)  in  the  p  metric.  To  this  end  we  prove  the  following  result. 

Theorem  5.1  Let  F  be  a  weakly  compact  subset  of  L2(Q),  with  Q  compact,  and  let  Vr(Q) 
be  the  family  of  probability  distribution  functions  on  Q  generated  by  F  as  a  set  of  densities 

VAQ)  =  {FeVmF'  =  f,feF}. 

Then  Vp(Q)is  a  p  compact  subset  of  ( V(Q),p )  where  p  is  the  Prohorov  metric  on  the  set 
V(Q)  of  all  probability  density  functions  on  Q. 

Proof:  Let  {Fh}  G  7 'V(Q),  where  Fn  —  f  fnds,  fn  G  F .  Since  weak  compactness  of  F 
implies  that  F  is  norm  bounded  we  have  that  {/„}  C  F  is  weakly  sequentially  compact,  i.e., 
there  exists  an  /  G  F  and  a  subsequence  {fnk}  of  {/„}  such  that  frik  — ^  /  with  ||/nJ|i,2  <  M. 
Then  for  any  continuous  function  g  G  C(Q)  we  have 


which  in  turn  implies  that  Fnk  — >  F  —  f  fds  in  the  p- metric.  Thus  Vp(Q)  is  sequentially 
compact.  Next  we  establish  that  this  set  is  also  closed  in  the  p-metric.  Let  {Fn}  C  Vt(Q) 
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be  such  that  Fn  — >  F  G  'P(Q)  in  the  p-metric.  We  must  argue  that  F  G  Vr(Q).  Since 
Fn  =  f  fnds,  we  have  that 

f  gfndx  ->  f  gdF  (5.2) 

for  every  p  G  C(Q).  Since  {/„}  C  T  we  can  assert  as  above  that  there  exists  {/nfc}  and 
/  G  T  such  that 

/  9fnkdx  -»•  /  p/cfe 

for  all  g  G  C(Q),  which  implies  that 


for  every  h  G  L2(Q),  since  is  bounded  in  L2(Q)  and  C'(Q)  is  dense  in  L2(Q).  Now  by 

(5.2) 

[  gdF  =  [  gfdx  for  all  g  G  L2(Q), 

JQ  JQ 

i.e. ,  F  is  absolutely  continuous  with  F'  =  /,  which  implies  that  F  G  T .  Thus  we  can  conclude 
that  Vf(Q)  is  compact  in  (V(Q),p). 


This  result,  combined  with  the  continuous  dependence  statement  above,  guarantees  that  the 
inverse  problem  is  solvable,  i.e.,  a  minimizer  of  the  cost  functional  J(P)  exists  over  the  set 
Vjr{Q)  C  V(Q).  We  note,  that  although  the  previous  computational  framework  utilizing 
Dirac  measures  is  valid  here,  it  may  be  desirable  to  develop  ‘smoother’  approximations  to 
elements  of  Vr{Q)-  Namely,  suppose  that  /  G  T  and  F  G  V(Q)  with  F  =  f  f.  Since 
T  C  L2(Q),  we  can  formulate  a  piecewise  linear  spline  approximation  to  /,  i.e.,  let 

N 

J— o 

(where  the  lj- s  are  the  usual  piecewise  linear  splines)  such  that  fN  — >  /  in  L2(Q),  and  the 
b^-s  are  rational  numbers.  This  implies  that 

f  dfNdx  -*•  f  gfdx 

Jq  Jq 

for  all  g  G  L2(Q),  hence  for  all  g  G  C(Q),  which  yields 

0. 

Let  FN  =  {h  G  L2(Q)\h(x)  =  Y,f=obf£f(x)},  where  the  b^-s  are  rational.  Thus  we  can 
conclude  that  the  set 

V(Q)  =  {F£P(Q)\F=  J  f,f£U?TN]  (5.3) 

is  dense  in  VriQ)  in  the  p-metric.  Hence  piecewise  linear  splines  provide  an  alternative 
way  to  approximate  elements  of  Vp{Q)  in  our  computational  work.  (Recall  that  elements 
of  U VM(Q),  which  is  a  dense  subset  of  V(Q),  can  also  be  used  to  approximate  elements  in 
Vt{Q)  even  though  the  elements  of  VM(Q)  are  not  in  Vt{Q)-) 
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6  A  computational  example 


Figure  2:  Comparison  of  the  optimized  model  normalized  acceleration  with  two  internal 
variables  (dashed  line)  and  Verburg  acceleration  data  (solid  line) 

In  this  section  we  present  an  inverse  problem  example  to  illustrate  the  computational 
aspects  of  the  theory  developed  above.  The  particular  example  presented  is  part  of  the  com¬ 
putational  effort  reported  in  [2].  “Data”  was  simulated  using  the  experimental  frequencies 
of  Verburg  [27]  for  a  medium  with  mechanical  properties  similar  to  the  lung  tissue  located 
between  the  heart  valves  and  the  chest  wall.  We  shall  refer  to  this  as  the  Verburg  data; 
details  on  precisely  how  this  data  is  generated  can  be  found  in  [2]  but  we  note  it  contains 
no  damping  or  dissipation.  It  was  used  as  observations  for  shear  acceleration  a*  at  in 
a  least  squares  formulation  similar  to  (5.1).  The  model  chosen  was  of  the  form  (2.1)-(2.3) 
with  Co  =  0  and  <re  given  by  (2.4)  with  7  =  —(3  =  —1.0  taken  fixed.  In  the  optimization 
problem,  the  class  of  measures  was  restricted  to  the  set 

V\Q)  =  [P  6  V(Q)\P(t)  =  -2Sn(T)  +  t)}, 

where  in  this  case  Q  was  taken  as  a  Cartesian  product  of  compact  intervals  of  positive  r 
values,  i.e. ,  Q  =  [rf ,  rf1-]  x  [r^,  In  the  example  we  also  estimated  Ci,  C2  in  (2.2)  as  well 
as  a  in  (2.4).  Thus  the  search  was  carried  out  over  five  Euclidean  parameters.  The  best  fit 
parameters  found  were  C{  =  44.78,  C 2  =  96.023,  t*  =  .000416,  r2*  =  .0038,  a*  =  .24.  The 
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model  simulation  with  these  optimal  parameters  is  compared  to  the  Verburg  data  in  Figure 
2  where  comparison  plots  are  given  in  both  the  time  domain  and  the  frequency  domain. 

We  note  that  one  could  also  use  elements  of  V{Q)  of  (5.3)  in  place  of  the  elements  from 
VM(Q)  in  these  problems.  While  such  an  approximation  would  yield  a  system  with  a  con¬ 
tinuous  spectrum  of  relaxation  times,  the  computational  efforts  in  this  class  of  problems 
promises  to  require  somewhat  more  substantial  efforts  than  use  of  the  elements  in  VM(Q). 
Since  the  elements  in  V2(Q)  performed  quite  well  with  the  Verburg  data,  we  did  not  pursue 
piecewise  linear  spline  approximations  for  these  problems.  The  authors  of  [2]  did  consider 
certain  parameterized  nonlinear  versions  of  (2.2),  i.e.,  the  y£j  term  is  replaced  by  a  parame¬ 
terized  nonlinearity  g(Sintlb )  consisting  of  piecewise  linear  splines  in  the  internal  strain,  but 
they  noted  no  improved  performance  over  that  when  using  the  VM {Q)  approximate  systems 
with  the  Verburg  data.  We  note  here  that  these  nonlinear  systems  considered  in  [2],  which 
can  be  written 

+g(£int',b)  =  C^ae(ux(t)), 

can  also  be  readily  treated  by  the  probabilistic  framework  developed  here.  In  this  case,  one 
can  take  measures  over  the  parameters  b  and  C,  where  £int(f)  =  £int(t-,b,C),  instead  of  over 
the  parameters  r  and  C. 

7  Concluding  remarks 

In  the  above  discussions  we  have  presented  a  theoretical  framework  based  on  a  probabilistic 
multiscale  formulation  for  molecular  based  hysteresis  in  tissue-like  materials.  These  efforts 
were  motivated  by  and  formulated  in  the  context  of  shear  propagation  in  which  reptation  for 
long  strand  molecules  plays  a  fundamental  role  and  where  the  molecular  dynamic  parameters 
vary  across  populations  of  molecules  in  an  unknown,  probabilistic  manner.  We  developed 
an  inverse  problem  approach  for  measure  dependent  dynamical  systems  (in  this  case,  partial 
differential  equations  for  shear  displacements)  by  making  use  of  the  Prohorov  metric  topology 
on  measures.  The  formulation  offers  a  computationally  tractable  alternative  to  the  well- 
known  Fung  kernel  approach  for  viscoelastic  biomaterials. 

While  the  motivating  applications  are  important  in  their  own  right,  we  are  confident 
that  the  underlying  ideas  have  potential  for  a  much  wider  class  of  applications.  In  addition 
to  the  areas  mentioned  in  the  Introduction,  we  feel  that  such  an  approach  will  eventu¬ 
ally  prove  relevant  in  diverse  areas  such  as  electromagnetic  polarization  and  conductivity 
in  heterogeneous  materials  (see  [9], [12]  and  the  references  cited  therein),  industrial  poly¬ 
meric  melts  ([16], [17], [18], [24]),  as  well  as  in  other  biological  applications  [4]  (e.g.,  disease 
pathogenesis,  epidemiology,  ecological  migrations,  and  genomic  to  system  response  models- 
bioinformatics!).  The  fact  that  the  theoretical  formulation  leads  readily  to  both  computa¬ 
tionally  and  theoretically  sound  approximations  significantly  enhances  the  attractiveness  of 
this  approach. 
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